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The  Huygen-Fresnel  formulation  for  analysis  of  the 
propagation  of  an  electromagnetic  wavefront  is  extended  from 
the  monochromatic,  steady  state  condition  to  the  transient 
condition.  The  propagation  of  a  plane  wave  and  the  transient 
response  of  a  simple  lens  are  analyzed  to  illustrate  the 
application  of  the  technique  and  to  show  that  the  results  are 
consistent  with  those  expected  from  the  more  familiar  steady 
state  formulations. 


TRANSIENT  RESPONSE  AND  HUYGEN-FRESNEL 


Introduction 

Calculations  of  diffraction  patterns  and  analysis  of  other 
propagation  problems  by  Huygen-Fresnel  wavefront  construction  techniques 
are  common  content  of  optical  texts.1  It  is  generally  recognized  that 
the  results  achieved  are  approximate,  but,  for  many  sets  of  conditions, 
sufficiently  accurate  to  make  the  techniques  a  valuable  engineering  tool. 
In  the  present  paper  the  basic  Huygen-Fresnel  formulation  is  extended 
from  the  monochromatic,  steady  state  condition  to  the  transient  condition. 
The  steady  state  condition  is  a  special  case  derivable  from  the  transient 
results.  My  motivation  for  offering  this  development  is  based  upon  my 
belief  that  in  some  cases  the  approach  yields  increased  insight  into 
expected  results,  and  offers  mathematical  simplification  as  compared  to 
the  steady  state  formulation  under  some  conditions.  The  approach  also 
yields  insight  into  the  nature  and  magnitude  of  error  which  may  be 
attributable  to  the  approximations  inherent  in  the  Huygen-Fresnel 
formulation. 

Linearity  and  Superposition 

The  development  in  this  paper  is  based  upon  concepts  of  linearity 
and  the  resulting  superposition  theorems.  These  fundamental  notions 
have  been  around  for  a  very  long  time2,  and  are  briefly  reviewed  to  set 
the  stage  for  the  developments  in  this  paper.  The  concepts  are  sum¬ 
marized  in  Fig.  1.  Line  (a)  depicts  a  general  input  to  the  medium  f(t) 


A 


MEDIUM 


(a) 


F(t) 


(b) 


G(t) 


(c) 


F(t)  +  G(t) 


(d) 


A(&j)sin[<ut  +  qSloi)] 


s(t) 

(e)  - 


Fig.  1.  Concepts  of  Linearity  and  Superposition 
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which  produces  an  output  F(t).  It  is  assumed  that  the  medium  is  in¬ 
variant  in  time  over  the  intervals  of  Interest.  This  is  equivalent  to 
stating  that  the  same  input  function  inserted  at  a  different  time, 
f(t-x),  will  produce  an  identical  output  simply  delayed  in  time  F(t-r). 
Line  (b)  depicts  another  general  input,  g(t),  which  produces  an  output, 
G(t).  The  critical  test  of  linearity  lies  in  the  superposition  indicated 
in  line  (c),  i.e.,  for  a  linear  system  an  input  consisting  of  the  sum 
of  two  or  more  functions  will  produce  an  output  which  is  the  sum  of  the 
outputs  which  would  be  produced  by  each  of  the  inputs  individually. 

It  is  this  linearity  of  superposition  which  makes  it  possible  to  disect 
the  input  into  arbitrarily  chosen  components,  as  for  example,  individual 
sinusoidal  components  obtained  by  Fourier  decomposition  of  the  input  as 
depicted  in  line  (d) .  Since  any  input  can  be  synthesized  by  Fourier 
analysis,  the  properties  of  the  linear  system  can  be  specified  by  noting 
the  output  amplitude  and  phase  shift  for  a  unit  amplitude,  zero  phase 
input  for  each  of  all  possible  frequencies  in  the  entire  spectrum.  Since 
a  Dirac  function  6 ( t )  as  shown  on  line  (e)  contains  all  possible 
frequencies  with  fixed  phase  relationship,  the  output,  s(t),  from  such 
an  input  can  be  Fourier  analyzed  to  determine  the  transfer  function, 
i.e.,  the  attenuation  and  phase  shift  associated  with  each  frequency 
throughout  the  spectrum.  Since  any  input  function  can  also  be  decom¬ 
posed  into  an  infinite  set  of  Dirac  functions  of  appropriate  amplitude 
and  time  position,  the  output  from  any  generalized  input  will  be  the 
convolution  of  the  input  with  the  transient  response,  s(t). 


Propagation  of  a  Plane  Wave 


The  propagation  of  a  plane  wave  is  a  simple  case  to  treat  and  one 
for  which  we  have  strong  intuitions  as  to  the  proper  answer.  Consider 
the  following  geometry: 


y 


In  the  x-y  plane  there  is  a  plane  monochromatic  wavefront.  The  first 
step  will  be  to  use  a  Huygen-Fresnel  calculation  to  determine  the  ampli¬ 
tude  and  phase  of  the  field  at  point,  P. 

The  increment  of  field  at  P  due  to  the  incremental  area  dS  is, 
dU(P)  -  -  j  Q(9)  ei2W  X  e"l[wt+4,(a>)]  dS.  (1) 

The  term, 

K<0)  -  -  j  Q(6)  ,  (2) 

is  the  inclination  factor3,  \  is  the  wavelength,  and  w  is  the  angular 
frequency.  A(<d)  is  the  amplitude  of  the  monochromatic  wave  in  the  x-y 
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plane  and  $(w)  is  its  phase  in  that  plane.  By  noting  that. 


where  C  is  the  velocity  of  light,  eq.(l)  becomes. 


dU(P) 


i  u> 

2irC 


Q(0> 


A(o>)  iw  c  -i[wt-K}>(w)].  dS 


(3) 


(4) 


Therefore, 


vm  .J  2i|l  |.lu  *<„)  *<“j  ds  . 


(5) 


Equation  (5)  gives  the  amplitude  and  phase  of  the  monochromatic 
wave  at  point  P  generated  by  a  monochromatic  plane  wave  in  x-y  having 
an  amplitude  A(to)  and  a  phase  of  Assume  now  that  the  electro¬ 

magnetic  disturbance  in  x-y  is  not  a  monochromatic  wave  but  instead  is 
a  temporal  signature,  G(t) .  In  terms  of  monochromatic  components. 


G(t) 


+  .°° 


J  A(u>) 


e~i  [tat  +  <Kw)] 


dw 


(6) 


To  obtain  the  temporal  signature  at  point  P, 
+  “ 

V(t)  -  |  U(P)  du  . 

_  00 

Substituting  eq.(5)  into  eq.(7)  gives, 


(7) 


v<t)  -j  j  /-  iu  am  '  t  >  +  *(“)!J  as*.  • 

-  oo  fi  ^ 


(8) 
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Interchanging  the  order  of  integration  gives. 


+  00  / 
/«£  I  -f 

fi  -  00  L 


la>  A(o)>  e 


-i[w(t  -  -r 


-  f  +  <Kw)]j 


doidS  .  (9) 


The  inner  integral  can  be  seen  to  be  of  the  form  of  eq.(6)  except  for  the 

g 

presence  of  the  -i^  and  the  time  delay  —  .  It  can  be  noted  by  inspection 
that, 

Dt  G(t  -  -  J  -  iw  A(u>)  e"i[!l){t  "  +  ^“^dw,  (10) 


where  Dt  is  the  derivative  with  respect  to  time.  Substituting  in  eq.(9) 


gives 


V(t)  *  J  2^  °t  G(t  ‘  dS 


Reversing  the  order  of  the  derivative  and  integral  yields. 


?(t>  '  »t  (^2^  °(t  '  !>  “S|  ' 


The  incremental  area  dS  will  be  chosen  to  be  a  ring  of  radius  W,  and 
thickness  dW.  Therefore, 


dS  -  2tt  W  dW  . 


t«  .  1  -  (.*2+r2)h 

c  c 
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Then, 


,  ,  w  dW  , 
dt'  =  -r- rr-rsH 
c(w2+r2) 


Combining  (15)  in  (13)  gives, 

dS  =  27TC  (w2+r2)^  dt'  , 

or. 


dS  =  2ttc  s  dt* 


(15) 


(16) 


(17) 


Substituting  this  result  in  eq.(12)  yields, 

00 

V(t)  =  D(t)  |  Q(6)  G(t-t')  dt'  .  (18) 

r/c 

To  gain  insight,  let  us  start  with  the  simple  but  invalid  assumption  that, 

Q(0)  5  1  .  (19) 

Let  me  also  define  a  function 

F(t')  =  0  t’  <  -  (20) 

c 

and 

F(t')  =  1  t'  >  -  (21) 

—  c 

so  that,  + 

V(t)  =  D(t)  f  F(t')  G(t-t’)  dt’  .  (22) 

_  00 

Equation  (22)  is  a  convolution  integral  and  indicates  that  to  obtain  V(t) 
we  first  convolve  the  input  function  of  time  G(t)  with  the  step  function 
F(t)  and  then  differentiate  the  resultant  function.  For  any  function  G(t) 
this  produces. 
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V(t)  =  G  (t  -  |)  .  (23) 

This  result  states  that  a  plane  wavefront  in  x-y  with  any  temporal  signature 
will  be  reproduced  at  point  P  with  no  alteration  in  either  the  amplitude  or 
the  temporal  signature  other  than  the  expected  delay  in  arrival  due  to  the 
finite  propagation  velocity.  This  result,  although  based  upon  the  invalid 
assumption  as  to  Q(0) ,  produces  an  answer  which  is  in  keeping  with  the 
intuitive  notion  as  to  how  plane  waves  are  propagated  in  a  non -dispersive 
medium. 

How  would  this  result  be  altered  by  making  a  "valid"  assumption  as  to 
to  Q(0)?  From  reference  3, 

K(0)  =  -  ~  (1  +  cos0)  .  (24) 

From  eq.(2)  and  (24), 

Q (0)  =  h  +  h  cos 0  .  (25) 


From  the  earlier  sketch, 


Let, 


and  from  (14) , 


cos©  =  — 
s 


t’ 


£ 

c 


then. 


COS0  =  — 
t' 


(26) 


(27) 


(28) 


(29) 
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so  that. 


Q(6)  “  ^  +  2t» 


By  analogy  to  eqs.  (20)  and  21)  I  will  redefine, 
F(t')  =  0, 


(30) 


t’<  - 
c 


and 


F(t’)  =4+2?, 


f 

c 


(31) 


(32) 


Returning  to  eq.(22)  with  the  redefined  F(t')  let  us  examine  the  special 
case  in  which, 

G(t-t')  *=  «(t-t’)  ,  (33) 

where  <5(t-t')  is  a  Dirac  function  having  value  only  at  t  -  t'.  Since, 


+  oo 


|  F(t')  6(t-t ' )  dt’  =  F(t)  , 


(34) 


then. 


Therefore 


V(t)  «  D(t)  F(t) 


V(t)  =  6(to)  -  ^ 


A  sketch  of  V(t)  is  as  follows: 


V(t) 


t  >  t 


(35) 


(36) 


t  •> 
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This  result  Indicates  that  a  Dirac  function  plane  wavefront  will  not 
remain  a  Dirac  function  after  propagation,  i.e.,  there  is  a  transient 
response  associated  with  the  propagation  process.  Note  that  the  Dirac 
function  portion  of  V(t)  is  the  result  achieved  when  Q(0)  was  assumed 
to  be  unity. 

The  Fourier  transform  of  V(t)  gives  the  steady  state  condition 
in  terms  of  the  amplitude  and  phase  for  all  possible  monochromatic  plane 
waves.  The  exact  form  of  this  Fourier  transform  is  unimportant  to  the 
present  arguments.  Certain  assymptotic  properties  of  the  transform  are 
important  and  can  be  deduced  easily.  To  do  so,  consider  first  the 
Dirac  function  portion  of  V(t).  The  Fourier  transform  of  the  Dirac 
function  is  a  flat  spectrum  of  unity  amplitude  with  a  phase  shift, 

ct 

♦  =  2n  -jS.  =  u  to  ,  (37) 

where  w  is  the  angular  frequency.  Consider  now  the  portion  of  v(t) 
involving  the  inverse  square  of  t'.  At  u>  »  o  the  contribution  to  the 
spectrum  from  this  component  will  be, 

oo  t 

B(u>  =  o)  =  -  J  dt  ,  (38) 

t 

o 

or, 

t  00 

B(uj  =  o)  =  2t  I  t  *  -  %  *  (39) 

o 

At  very  high  frequencies  this  portion  of  V(t)  will  have  a  spectrum 
which  approaches  zero.  The  combined  spectrum  is  therefore  as  shown 
labeled  (a)  in  the  sketch, 
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(b) 


whereas  the  line  labeled  (b)  is  the  result  when  Q(0)  is  assumed  to  be 
unity . 

As  a  practical  matter,  it  is  Important  to  know  the  conditions 
under  which  (a)  and  (b)  become  approximately  equal.  This  will  occur 
when, 

t  »  —  .  (AO) 

o  to 

This  is  equivalent  to  stating  that. 


or. 


r  »  X  . 


(41) 


(42) 


Therefore  when  r  »  X  the  result  will  be  numerically  similar  for  the 
"valid"  and  "invalid"  choices  of  Q(8). 
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or 


T(r)  -  -  WR!+i  2  +  Vl2+r  2  -  V*t2+r2  -  Vl2-r2 


(44) 


For  simplicity  assume  that  I»r  and  R  >>rm  80  that  by  binomial  expansion. 


x(D  -  f 


R+ 


2  2  ] 

2R  +  1  +  21  R  2R  1  21  J  ’ 


(45) 


or. 


T(r)  - 


^1  „  ifl  .  l)  „2  (l  .  l) 

F  r®  (r  +  1/  ~  r  U  +  I> 


(46) 


or 


T(t)  -fe(rmZ-t!)(5  +  l) 


(47) 


Since, 


I  +  I-I 

R  I  f 


(48) 


where  f  is  the  lens  focal  length,  then, 
x(r) 


2c?  (rm2-lZ) 


(49) 


In  contrast  to  the  case  of  the  propagation  of  a  plane  wave,  the 
properties  of  the  simple  lens  are  such  that  the  on-axis  imaging  of  a  point 
source  Dirac  function  the  F  function  is, 


F(t’)  -  S(tQ) 

and 

F(t’)  -  0  ,  t  1»  to 


(50) 


(51) 
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When  convolved  with  the  input  time  function  here  assumed  to  be  a  Dirac 
function  and  differentiated  in  the  manner  of  eq.(22),  the  resultant 
transient  response  is  a  differentiated  Dirac  function,  i.e.,  the  derivative 
operator.  In  other  words,  a  Dirac  function  source  produces  an  image  which 
is  not  a  Dirac  function  but  the  derivative  of  a  Dirac  function. 

The  power  spectrum  of  this  on  axis  image  is  found  by  Fourier 
transforming  the  transient  response  and  squaring  the  amplitude  at  each 
frequency.  The  derivative  operator  has  a  Fourier  transform  which  is  of 
the  form, 

I A ( f ) |  =  Kf,  (52) 

where  f  is  frequency  and  K  is  a  constant.  The  power  spectrum  is  therefore 
of  the  form, 

I A(f ) |  -  K2f2  .  (53) 

In  terms  of  wavelength, 

c  ^ 

I A(f ) | 2  -  K  J2  .  (54) 

This  displays  the  inverse  square  of  wavelength  which  we  recognize  as 
characteristic  of  the  on-axis  imagery  of  a  point  source  by  a  diffraction 
limited  optical  system. 

Consider  now  a  point  in  the  image  plane  off-axis  in  the  X  direction 
by  a  distance  AX.  The  distance  between  that  point  and  any  point  x',y'  in 
the  lens  plane  is 

P  -  V  (x'-Ax)2  +  y’2  +I2  . 


(55) 
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If  I»  (x'-Ax)  and  y'  then. 


P  .  z  +  (s'-*)2^'2 

r  i  2J 


(56) 


The  total  transit  time  from  an  on-axls  Dirac  function  point  source  Is 
therefore. 


T  .  I  /  „  +  fj-  +  I  + 


\  8  +  it  +  i  +  - 1 


(57) 


Expanding,  collecting  terms  and  simplifying  gives. 


T  -  -  (R  +  I  +  — 


rm\,  1  ItA  1  J"x,Ax\ 

*/  ■I’Tni  1  /  * 


(58) 


The  equation  has  been  arranged  into  three  terms  because  each  has  special 
significance.  The  first  term  is  the  transit  time  from  the  point  source, 
through  the  center  of  the  lens  to  the  center  of  the  image  plane.  The 
second,  is  a  term  which  accounts  for  the  quadratic  phase  shift  in  the 
image  plane  familiar  to  conventional  optics.  The  third  term  is  the  only 
one  containing  the  spatial  variables  of  the  lens,  i.e.,  x'.  The  first  two 
terms,  therefore,  represent  only  time  delays  because  they  are  independent 
of  lens  plane  coordinates.  The  third  term  is  a  transient  response  because 
each  x'  position  on  the  lens  has  a  different  arrival  time  at  the  image 
plane  location  Ax,o.  Specifically  in  a  time  interval  dx  there  will  be  a 
return  from  a  region  of  lens  such  that, 


dx 


Ax 

Cl 


dx' 


(59) 
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The  signal  will  be  received  in  the  image  plane  at  Ax,o  within  a  time  interval. 


(60) 


f 

r  2 

R  +  T  +  — 

\ 

\  +  I  < 

M 

1,1  ll 

R  +  1  +  2f 

21  1 

[  C  \ 

V.  i 

1  1 

L  j 

1  V, 

Since  dx'  is  independent  of  x',  the  amplitude  of  the  return  will  be  determined 
by  the  height  of  the  lens  at  each  x'  position.  For  the  circular  lens. 


■Yr  2-x’2 
m 


(61) 


The  amplitude  of  the  signal  is  therefore  proportional  to. 


y’dx'  -  Vr  2-x'2  -r^-  dT 
m  Ax . 


(62) 


For  simplicity  let  the  sum  of  the  first  two  terms  of  eq.(58)  be  T  ,  then. 


T  -  T 


■‘W  • 


or. 


x'  =  -T~  {t  -t)  . 

OX  0 


Substitution  in  eq.(62)  gives, 
y’dx’  =  /r  2  - 


/T  \  (t  -t}2  *T“  dT 
m  ( Ax ;  o  Ax 


From  eq.(60),  the  range  of  x  is, 


o  c  ^  I  j 


(63) 


(64) 


(65) 


(66) 
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A  sketch  of  eq.(65)  is  as  follows: 


r  Cl 


The  Fourier  transform  of  this  function  Is  the  modified  Bessel  function,  i.e.. 


G  =* 


Kirr  2 
m 


(67) 


where  f  is  the  temporal  frequency,  and  K  is  a  constant  of  proportionality. 

As  described  earlier  in  this  writing,  the  derivative  of  eq.(25)  is  the 
transient  response.  A  derivative  in  the  time  domain  is  equivalent  to  multi¬ 
plication  by  frequency  coupled  with  a  90°  phase  shift.  Therefore, 


G'  -  iK  f  r  2 
l  m 


2frAxr  f 
_ m_ 

IC 


(68) 


The  following  observations  may  be  made  with  respect  to  eq.(68).  For  a  fixed 
f,  i.e.,  a  fixed  wavelength,  the  amplitude  as  a  function  of  Ax,  will  be  a 
modified  Bessel  function  therefore  predicting  the  classic  monochromatic 
diffraction  image  of  a  point  source  by  a  circular  lens.  It  can  also  be  noted 
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that  the  first  null  of  the  modified  Bessel  function  occurs  when  the  argument 
is  equal  to  1.22ir  so  that, 

Axr  f 

2it  — jS-  -  1.22ft  (69) 

For  Ax«I,  the  angular  position  of  Ax  in  the  image  plane  is. 


a 


Ax 
I  * 


and  since. 


(70) 


(71) 


then, 


or. 


a  r 


m 


=  0.61 


a  *  0.61 


m 


(72) 


(73) 


which  is  the  familiar  resolving  power  limitation  associated  with  a  diffraction 
limited  circular  aperture. 
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